#
#The Following R-Code creates a set of boxplots which demonstrate the effects of
#contiguity, alliance, joint democracy, and major power status on the probability
#of dyadic non-harmony for our ZiOPC replications of Senese (1997).
#
#These correspond to Figure 3
#


# clear memory
rm( list=ls() ) 

## set the memory to the max
memory.limit(size=4095)
memory.size(max = TRUE)

						     
# load some libraries			     
library(mvtnorm)
library(car)                                                          
library(scatterplot3d)
library(plotrix)
library(Hmisc)
library(Zelig)
library(foreign, pos=4)
library(pscl)

# set the seed
set.seed(2)

#Input the coefficient estimates from our Senese ZiOPC model
estimate<-cbind(0.74757, -3.77804, -1.97628, -0.20074, -2.76473,  3.59316,  0.29423, -0.15366,  1.04496, -0.78682, -0.20957, -0.44798,  0.09295,  0.06872,  0.10578)

#Input the variance-Covariance matrix from our Senese ZiOPC model
vcv<-matrix(0,nrow=15,ncol=15)
vcv[1,1:15]<-cbind(0.0883577000,0.0001827000,0.0001265000,-0.0010850000,-0.0087820000,-0.0510413000,0.0048020000,0.0007346000,0.0148600000,0.0849983000,0.0009193000,-0.0018640000,0.0015300000,0.0014770000,0.0366544000)
vcv[2,1:15]<-cbind(0.0001827000,0.0156100000,0.0008518000,0.0007242000,-0.0014870000,-0.0068620000,-0.0000898500,-0.0002846000,-0.0000894000,0.0007570000,-0.0002378000,-0.0001220000,0.0000574300,0.0002134000,0.0014110000)
vcv[3,1:15]<-cbind(0.0001265000,0.0008518000,0.0034400000,0.0006606000,-0.0013900000,-0.0063930000,-0.0000840800,-0.0002682000,-0.0000923700,0.0006555000,-0.0002276000,-0.0001109000,0.0000650500,0.0002120000,0.0013000000)
vcv[4,1:15]<-cbind(-0.0010850000,0.0007242000,0.0006606000,0.0014520000,-0.0010950000,-0.0037740000,-0.0001641000,-0.0002608000,-0.0002834000,-0.0007314000,-0.0001668000,-0.0000657800,0.0000586000,0.0001751000,0.0006350000)
vcv[5,1:15]<-cbind(-0.0087820000,-0.0014870000,-0.0013900000,-0.0010950000,0.0034410000,0.0105900000,-0.0004399000,0.0002836000,-0.0016370000,-0.0089090000,0.0001432000,0.0003950000,-0.0000588700,-0.0003107000,-0.0057300000)
vcv[6,1:15]<-cbind(-0.0510413000,-0.0068620000,-0.0063930000,-0.0037740000,0.0105900000,0.1389632000,-0.0047780000,-0.0016590000,-0.0078130000,-0.0596902000,0.0039100000,0.0033020000,-0.0009676000,-0.0023230000,-0.0263188000)
vcv[7,1:15]<-cbind(0.0048020000,-0.0000898500,-0.0000840800,-0.0001641000,-0.0004399000,-0.0047780000,0.0036130000,-0.0000570500,0.0008206000,0.0051350000,-0.0009321000,-0.0002150000,0.0000090420,-0.0000065210,0.0018120000)
vcv[8,1:15]<-cbind(0.0007346000,-0.0002846000,-0.0002682000,-0.0002608000,0.0002836000,	-0.0016590000,-0.0000570500,0.0054280000,0.0001020000,0.0010540000,-0.0000541400,-0.0025480000,0.0000272700,-0.0000181100,-0.0001715000)
vcv[9,1:15]<-cbind(0.0148600000,-0.0000894000,-0.0000923700,-0.0002834000,-0.0016370000,-0.0078130000,0.0008206000,0.0001020000,0.0032570000,0.0144100000,0.0001673000,-0.0003607000,0.0000208800,0.0000945300,0.0060810000)
vcv[10,1:15]<-cbind(0.0849983000,0.0007570000,0.0006555000,-0.0007314000,-0.0089090000,	-0.0596902000,0.0051350000,0.0010540000,0.0144100000,0.0842974000,0.0001716000,-0.0024040000,0.0001517000,0.0001292000,0.0360801000)
vcv[11,1:15]<-cbind(0.0009193000,-0.0002378000,-0.0002276000,-0.0001668000,0.0001432000,0.0039100000,-0.0009321000,-0.0000541400,0.0001673000,0.0001716000,0.0014360000,0.0000086240,-0.0000784400,-0.0001115000,0.0002176000)
vcv[12,1:15]<-cbind(-0.0018640000,-0.0001220000,-0.0001109000,-0.0000657800,0.0003950000,0.0033020000,-0.0002150000,-0.0025480000,-0.0003607000,-0.0024040000,0.0000086240,0.0037990000,0.0000763100,0.0001583000,-0.0009769000)
vcv[13,1:15]<-cbind(0.0015300000,0.0000574300,0.0000650500,0.0000586000,-0.0000588700,-0.0009676000,0.0000090420,0.0000272700,0.0000208800,0.0001517000,-0.0000784400,0.0000763100,0.0019740000,0.0015060000,0.0000399000)
vcv[14,1:15]<-cbind(0.0014770000,0.0002134000,0.0002120000,0.0001751000,-0.0003107000,-0.0023230000,-0.0000065210,-0.0000181100,0.0000945300,0.0001292000,-0.0001115000,0.0001583000,0.0015060000,0.0019500000,0.0002329000)
vcv[15,1:15]<-cbind(0.0366544000,0.0014110000,0.0013000000,0.0006350000,-0.0057300000,-0.0263188000,0.0018120000,-0.0001715000,0.0060810000,0.0360801000,0.0002176000,-0.0009769000,0.0000399000,0.0002329000,0.0173118000)



#calculate the probability of non-harmony for 0->1 change in contiguity
nsims<-10000
contig<-matrix(1,nrow=nsims,ncol=1)

  for(j in 1:nsims){
  bsim<-rmvnorm(1, mean=estimate, sigma=vcv)
  XB<-(1*bsim[,5]+0*bsim[,6]+0*bsim[,7]+0*bsim[,8]+0*bsim[,9])
  ZG1<-pnorm(XB)
  
  XB<-(1*bsim[,5]+1*bsim[,6]+0*bsim[,7]+0*bsim[,8]+0*bsim[,9])
  ZG2<-pnorm(XB)
  
  pzero<-ZG2-ZG1
  contig[j,1]<-pzero
}



#calculate the probability of non-harmony for 0->1 change in alliance
nsims<-10000
alliance<-matrix(2,nrow=nsims,ncol=1)

  for(j in 1:nsims){
  bsim<-rmvnorm(1, mean=estimate, sigma=vcv)
  XB<-(1*bsim[,5]+0*bsim[,6]+0*bsim[,7]+0*bsim[,8]+0*bsim[,9])
  ZG1<-pnorm(XB)
  
  XB<-(1*bsim[,5]+0*bsim[,6]+1*bsim[,7]+0*bsim[,8]+0*bsim[,9])
  ZG2<-pnorm(XB)
  
  pzero<-ZG2-ZG1
  alliance[j,1]<-pzero
}



#calculate the probability of non-harmony for 0->1 change in democracy
nsims<-10000
democracy<-matrix(3,nrow=nsims,ncol=1)

  for(j in 1:nsims){
  bsim<-rmvnorm(1, mean=estimate, sigma=vcv)
  XB<-(1*bsim[,5]+0*bsim[,6]+0*bsim[,7]+0*bsim[,8]+0*bsim[,9])
  ZG1<-pnorm(XB)
  
  XB<-(1*bsim[,5]+0*bsim[,6]+0*bsim[,7]+1*bsim[,8]+0*bsim[,9])
  ZG2<-pnorm(XB)
  
  pzero<-ZG2-ZG1
  democracy[j,1]<-pzero
}


#calculate the probability of non-harmony for 0->1 change in major power
nsims<-10000
majorpower<-matrix(4,nrow=nsims,ncol=1)

  for(j in 1:nsims){
  bsim<-rmvnorm(1, mean=estimate, sigma=vcv)
  XB<-(1*bsim[,5]+0*bsim[,6]+0*bsim[,7]+0*bsim[,8]+0*bsim[,9])
  ZG1<-pnorm(XB)
  
  XB<-(1*bsim[,5]+0*bsim[,6]+0*bsim[,7]+0*bsim[,8]+1*bsim[,9])
  ZG2<-pnorm(XB)
  
  pzero<-ZG2-ZG1
  majorpower[j,1]<-pzero
}


#graph the above results via-boxplots:

#set up the graphing space
par(mfcol=c(1,4))
par(cex.lab=1.25)
par(cex.axis=1.25)
par(cex.main=1.25)

#begin the boxplots
boxplot(contig, xlab="Contiguity",ylab="Change in Pr(S=1)", main="Model: ZiOPC")
par(new=FALSE)
boxplot(alliance, xlab="Alliance",ylab="Change in Pr(S=1)",main="Model: ZiOPC")
par(new=FALSE)
boxplot(democracy, xlab="Joint Democracy",ylab="Change in Pr(S=1)",main="Model: ZiOPC")
par(new=FALSE)
boxplot(majorpower, xlab="Major Power",ylab="Change in Pr(S=1)",main="Model: ZiOPC")


